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We perform a von Neumann stability analysis on a common discretization of the Ein- 
stein equations. The analysis is performed on two formulations of the Einstein equations, 
namely, the standard ADM formulation and the conformal-traceless (CT) formulation. 
' The eigenvalues of the amplification matrix are computed for flat space as well as for 

. a highly nonlinear plane wave exact solution. We find that for the fiat space initial 

' data, the condition for stability is simply < 1- However, a von Neumann analysis 

^N) ^ for highly nonlinear plane wave initial data shows that the standard ADM formulation 

is unconditionally unstable, while the conformal-traceless (CT) formulation is stable for 
^\ 0.25<#i<l. 
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> ■ I. INTRODUCTION 

i> : 

^ . 

, Numerically solving the full 3D nonlinear Einstein equations is, for several reasons, a daunting task. 

QQ ■ Still, numerical relativity remains the best method for studying astrophysically interesting regions of the 

\ solution space of the Einstein equations in sufficient detail and accuracy in order to be used to interpret 

. measurements made by the up and coming gravitational wave detectors. Even though numerical relativity 

' is almost 35 years old, some of the same problems faced by researchers three decades ago are present today. 

O . Aside from the computational complexity of implementing a numerical solver for the nonlinear Einstein 

O^' equations, there exist several unsolved problems, including the well-posedness of certain initial value 

^ \ formulations of the Einstein equations and the proper choice of gauge. Not the least of these problems is 

bX)' numerical stability. A common thread in numerical relativity research over the past three decades is the 

^ I observation of high frequency (Nyquist frequency) noise growing and dominating the numerical solution. 

Traditionally, numerical studies have been performed with the initial value formulation of the Einstein 
equations known as the ADM 3-1-1 formulation [Q, in which the 3-metric and extrinsic curvature are the 
dynamically evolved variables. 

Lately, a formulation based on variables in which the conformal factor of the 3-metric and the trace 
of the extrinsic curvature are factored out and evolved separately has been studied. This conformal- 
traceless (CT) formulation was first introduced by Nakamura and Shibata Q and later slightly modified 
by Baumgarte and Shapiro The stability properties of the CT formulation were shown in |^ to be 
better than those of the ADM formulation for linear waves. The improvement in numerical stability J_n 
the CT formulation versus the ADM formulation was demonstrated in strong field dynamical cases in 
A step toward understanding the improved stability properties of the CT formulation was taken in 
where it was shown by analytically linearizing the ADM and CT equations about flat space that the CT 
system effectively decouples the gauge modes and constraint violating modes. It was conjectured that 
giving the constraint violating modes nonzero propagation speed results in a stable evolution. 

Here, we take another step towards understanding the improved stability properties of the CT system 
by performing a von Neumann stability analysis on discretizations of both the ADM and CT systems. 
We are led to the von Neumann stability analysis by Lax's equivalence theorem j^], which states that 
given a well posed initial value problem and a discretization that is consistent with that initial value 
problem (i.e., the finite difference equations are faithful to the differential equations), then stability is 
equivalent to convergence. Here, the words "stability" and "convergence" are taken to mean very specific 
things. Convergence is taken to mean pointwise convergence of solutions of the finite difference equations 
to solutions of the differential equations. This is the piece de resistance of numerical relativity. After all, 
what we are interested in are solutions to the differential equations. Stability, on the other hand, has a 
rather technical definition involving the uniform boundedness of the discrete Fourier transform of the finite 
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difference update operator (see |0] for details). In essence, stability is the statement that there should 
be a limit to the extent to which any component of an initial discrete function can be amplified during 
the numerical evolution procedure (note that stability is a statement concerning the finite difference 
equations, not the differential equations). Fortunately, the technical definition of stability can be shown 
to be equivalent to the von Neumann stability condition, which will be described in detail in the next 
section. 

While one cannot apply Lax's equivalence theorem directly in numerical relativity (the initial value 
problem well-posedness assumption is not valid for the Einstein field equations in that the evolution 
operator is not, in general, uniformly bounded ), numerical relativists often use it as a "road map" ; clearly 
consistency and stability are important parts of any discretization of the Einstein equations (curiously, 
convergence is usually implicitly assumed in most numerical studies). Code tests, if done at all, usually 
center around verifying the consistency of the finite difference equations to the differential equations (as 
an example of the extents to which some numerical relativists will go to check the consistency of the 
finite difference equations to the differential equations, see, e.g., §). Stability, on the other hand, is 
usually assessed postmortem. If the code crashes immediately after a sharp rise in Nyquist frequency 
noise and/or if the code crashes sooner in coordinate time at higher resolutions, the code is deemed 
unstable. 

We suggest that the stability of the code can be assessed before (and perhaps even more importantly, 
while) numerical evolutions take place. As will be seen in the next section, the stability properties of 
any given nonlinear finite difference update operator depend not only on the Courant factor but 
also on the values of the discrete evolution variables themselves. Therefore, during numerical evolutions 
of nonlinear problems, as the evolved variables change from discrete timestep to timestep, the stability 
properties of the finite difference operator change along with them! Ideally, one would want to verify 
that the finite difference update operator remains stable for each point in the computational domain at 
each timestep. While the computational expense of this verification would be prohibitive, verification at 
a reasonably sampled subset of discrete points could be feasible. 

The remainder of the paper is outlined as follows. Section || will describe the von Neumann stabil- 
ity analysis for a discretization of a general set of nonlinear partial differential equations in one spatial 
dimension. Included will be results from a von Neumann stability analysis of the linear wave equation 



discretized with an iterative Crank-Nicholson scheme. Section III will present the ADM and CT formula- 
tions of the Einstein equations restricted to diagonal metrics and further restricted to dependence on only 
one spatial variable. The von Neumann stability analysis is performed on iterative Crank-Nicholson dis- 



cretizations of both formulations with flat initial data. Section IV will repeat the von Neumann stability 



analysis from section [II with a nonlinear wave solution for initial data. 



II. THE VON NEUMANN STABILITY ANALYSIS 

Let a set of partial differential equations be given as 

^^Siu,d,u,d,'u) (1) 

where u = u{z,t) is a vector whose components consist of the dependent variables that are functions of 
the independent variables Q < z < N and t > (while we ignore boundary conditions in this paper, as 
we consider only interior points of equations with finite propagation speeds, the von Neumann analysis 
presented here is easily extended to include boundary treatments Q]). 
Now, consider a discretization of the independent variables z and t: 

=jAz, j = 0,l,2,...,m (2) 
tn^nAt, n = 0,1,2,... (3) 

along with a discretization of the dependent variables u(z, t): 

u'l ^u{z,,tn). (4) 

Furthermore, consider a consistent discretization of the set of differential equations in the form of Eq. |] 
that can be written in the form 
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Uj'^ —S^j{...,u'J_2,i!^_i,i!^,u^^i,u'j^2J---)- (5) 

As we will only be analyzing the iterative Crank-Nicholson discretization scheme, we assume a two-step 
method as written in Eq. |[ However, a von Neumann analysis does not depend on this, and could easily 
be performed for three-step methods such as the leapfrog scheme. 
Now, assume some initial data 

< j < m (6) 
for the discrete variables is given at time t ~ 0. The amplification matrix at Zj is given by 

rn o nn 
"^3 2^ dun 

The condition for numerical stability is that the spectral radius of the amplification matrix G" be 1 
or less for each wavenumber kz- That is, each eigenvalue of the amplification matrix G" should have 
a modulus of 1 or less for each discrete mode kz 

\\^\ < 1. (8) 

First, notice that for linear finite difference equations (Eq. ||), the amplification matrix G" does not 
depend on either the initial data it^ nor on the spatial index j. It only depends on the discretization 
parameters Az and At as well as the mode wavenumber kz- That is, for the linear case, once one specifies 
the discretization parameters Az and At, one only needs to verify that the eigenvalues of the amplification 
matrix G" have a modulus of less than or equal to 1 for all wavenumbers kz to insure the numerical 
stability of a numerical update operator, regardless of initial data Contrast this with the nonlinear 
case where the amplification matrix G" is not only a function of discretization parameters Az, At and 
the mode wavenumber kz, but also of initial data itj. In this case, not only must one verify that the 
amplification matrix G" have a spectral radius of 1 or less for each spatial index j (assuming the initial 
data i^j depends on j, which it will in general), but one must carry out this verification at every time step 
n, as the data iPj will, in general, change with increasing n. So, in the nonlinear case, the amplification 
matrix G" will depend on both the spatial index j and the temporal index n. In principle, one must 
verify that G" have a spectral radius of 1 or less for each mode wavenumber kz , for each spatial index j 
and for each time step n, in order to be confident that the solutions of the finite difference equations are 
converging to solutions of the differential equations (which is, after all, what we are ultimately interested 
in). 



A. Example: Linear Wave Equation 



In a von Neumann analysis of the advection equation was presented for the iterative Crank- 
Nicholson scheme. Here, as a prelude to computing a von Neumann analysis for discretizations of the 
equations of general relativity, we present a von Neumann analysis for a 2-iteration iterative Crank- 
Nicholson discretization of the wave equation in 1 dimension 



dt^ 
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We write the wave equation in first-order form as in Eq. 0. Defining D = H = and 



D 
H 

the wave equation in one spatial dimension becomes 

H 

dn 



(9) 



(10) 



du 
'dt 



an 



(11) 
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The iterative Crank-Nicholson discretization procedure begins by taking a FTCS (Forward Time Centered 
Space) step which is used to define the intermediate variables ^^■'m,^ : 



(1) 



— n+l 



At 



(n;Vi-n5'_i)/(2Az) 

(i?;Vi-^"-i)/(2Az) 



This intermediate state variable is averaged with the original state variable 



(l)^n+l/2 1/(1) =."+1 



(12) 



(13) 



which, in turn, is used to calculate '■^•'t/"'*'^, the state vector produced from the first iteration of the 



iterative Crank-Nicholson scheme: 



(1)^^+1 



(l)jjn+l/2 
3 



((i)n;^T-(^^n^T)/(2A.) 



(14) 



The averaged state 



(2),-7"+l/2 



is calculated 



and used to compute the final state variable '^^^u^^^ 



(2)^+1 



At 



(2)nj+i/2 

((2)n;+'/'-(4;+'/')/(2Az) 

((2)15";V2_ (2) ^"+1/2^/(2^^^ 



(15) 



(16) 



Using this second iteration as the final iteration, we can write Eq. ^ explicitly as the three equations 

2 A+ ^2 
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-(^;V2-^^,) + A^(l--)^7 



(17) 



(1 



n 



-n--i) 



(i-^)n;' + y(n: 



J+2 



n 



32 



(18) 



(19) 



where we have denoted a as the Courant factor At/Az. We compute the amplification matrix Eq. 0, 
which is independent of the spatial index j and the time index n: 



G = 



Gil Gi2 Gi3 
G21 G22 G23 

G31 G32 6*33 



(20) 



where the components are 



Gil 
G21 
G12 

Gi3 

G22 
G23 



1 

G31 







—3— sm t/fc 



Ml- V) 
G33 = (1 - 

G32 = -^sin(30fe)-z(T(l 



, g cos(26'fc) 
^) + ^cos(20,) 



(21) 
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where Ok = kz Az. The modulus of the eigenvahies of G are easily calculated and found to be (ignoring 
the eigenvalue that is identically 1) 



|Ai| = IA2I = y 1 - ^CT^sin^^fe + —aesm'^ek (22) 

By inspection, the von Neumann condition for stability, i.e. |A| < 1 for all kz, is a = < 2. However, 
it is instructive to look at the 9k dependence of the stability criterion. Figure |l| shows a plot of Eq. H] as 
a function of dk = kz Az for various values of the Courant factor a. As can be seen, the first eigenmode 
that goes unstable (i.e., has an eigenvalue whose modulus is greater than 1) as a is increased is the mode 
6k = kz Az — TT. This corresponds to modes of wavelength 2 Az, i.e. Nyquist frequency modes, which 
are precisely the modes that usually crop up and kill numerical relativity simulations. 




0.5 1 < ' < 1 < ' < 1 
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k 

FIG. 1. The largest eigenvalue (ignoring eigenvalues that are identically 1) of the amplification matrix as a 
function of 9k = fez Az for the 2-iteration iterative Crank-Nicholson discretization of the scalar wave equation. 
Various values of the Courant factor a = At/ Az are shown. 



III. GENERAL RELATIVITY 



Here, we present the analytic equations for the ADM and CT formulations that wc will use to study 
the stability properties of numerical methods in general relativity. We will discrctize the general relativity 
evolution equations using the same discretization used in the previous section for discretizing the scalar 
wave equation, namely, a 2-iteration iterative Crank-Nicholson scheme. Before analyzing the stability of 
the discretization for non-linear waves, we will perform a von Neumann analysis for the discretizations 
of both the ADM and CT equations about flat space. 



A. ADM Equations in 1-D (plane symmetry) 



The form of the metric that we will use to study stability properties of discretizations of the ADM 
form of the Einstein equations is given by 

ds^ = —a^ dt^ + gxx dx^ + gyy dy^ + gzzdz'^, (23) 

where the lapse a and metric functions g^x, gyy, and gzz are functions of the independent variables z 
and t. To put the evolution equations in first order form (Eq. we introduce the extrinsic curvature 
functions Kxx, Kyy, and Kzz, each of which are also function of the independent variables z and t. The 
evolution equations for this ADM system is given in the form of Eq. |l| as 
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^g^i 

= -2aKii 

at 
dK 

dt 

dot 

- = ^aK (24) 

where the index pair ii takes on the values {xx, yy, zz}, the index j is summed over the vahies {x, y, z}, K 
denotes the trace of the extrinsic curvature K = K^xl^xx + Kyy/ gyy + K^z/ Qzz, denotes the covariant 
derivative operator compatible with the 3-metric, and Ru denotes the three components of the 3-Ricci 
tensor, which are given explicitly as 



1 / dgxx \ ^ 1 dgxx dg. 



yy 



Ryy 



'^gxxgzz \ dz J 4:gyygzz dz dz 

1 dgxx dgzz 1 9^ 

dz dz 2gzz dz'^ 

2 

dgyy \ 1 dgxx dgyy 



2 



^zz 



dz J 4:gxxgzz dz dz 
1 dgyy dgzz 1 d'^gyy 



JZZ 



^ dz dz 2gz 



(25) 



dgxx \ , ^ f dg. 



Jyy 



Ixx"^ \ dz J Agyy"^ \ dz 
1 dgxx dgzz , 1 dgyy dgzz 



24. This local condition 



'^gxxgzz dz dz 4:gyygzz dz dz 

1 d'^gxx 1 d^gyy 

'2gxx dz'^ 2gyy dz"^ 

Throughout this paper, we use the so-called "1 + log" slicing condition in Eq 
on the lapse has been used successfully in several recent applications [|l0|,^|. Moreover, it is a local 
condition, and thus, the von Neumann analysis remains local (this would be in contrast with a global 
elliptic condition on the lapse, e.g. maximal slicing, in which the calculation of the sum in the definition 
of the amplification matrix in Eq. |7| would have nonzero global contributions) . 

The initial value formulation is completed by specifying initial data that satisfies the Hamiltonian and 
momentum constraints, given respectively by 

H = ^R + K^ - KjkK^'" = (26) 
Mz = V.K^z - VzK (27) 



B. Conformal-Traceless (CT) Equations in 1-D (plane symmetry) 

The form of the metric that we will use to study stability properties of the discretizations of the CT 
form of the Einstein equations, as defined in is given by 

ds^ ^ -a^ dt^ + e^'''{gxxdx^ + gyy dy'^ + (jzzdz'^), (28) 

where the lapse a, the conformal function 0, and the conformal metric components gxx, gyy, and gzz are 
functions of the independent variables z and t. The determinant of the conformal 3-metric is identically 1. 
Instead of evolving the extrinsic curvature components as in the ADM formalism, the extrinsic curvature 
is split into its trace (K) and traceless (An) components: 

Ku^e'^'^i^guK + Au). (29) 

In addition, the conformal connection function F^, defined by 
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= -Ofg^"" (30) 

is also treated as an evolved variable. There are therefore 10 evolution equations that are in the form of 
Eq. |l|, and given explicitly as 

-— — — aK 
dt 6 



'2a Ai. 



at 

^ = aAjkA^'' + laK^ - V.V^a 
i9i 3 

^ = a(e-4^i?,, - ^guA^'^A^'' + '^g^^K^ + KA,, - 2A,,A\) - (31) 
e-4'^(P,I?,a - )^g^{Dp^a) 

at 3(7^2 oz (7^^ az 

9a 



dt 



-aK 



where the indices of Aij are raised and lowered with the conformal metric gij , the index on the covariant 
derivative operator D with respect to the physical metric is raised and lowered with the physical metric 
9ij — ^'^'^9ijj 8,nd are the Christoffel symbols related to the conformal metric gij. Note that the 

Hamiltonian constraint has been substituted in for the 3-Ricci scalar in the equations for ^ and 

Also, the momentum constraint has been substituted in the equation for This corresponds to the 

"Mom" system from [Q, and the a — l,m = l, £,=0 system from Note that the Ricci components 
Rii can be written in terms of the conformal Ricci components Ru: 

R^^ = R^^ + 2^^? - 2~g,,{V + ) + (32) 

dz g^^ dz 

''^z 4(— -2—), 33 
dz dz^ 

where the indices of the covariant derivative operator V with respect to the conformal metric are raised 
and lowered with the conformal metric. The conformal Ricci components can in turn be written as 

p 1 ^2.9, , dV' 

+ ^f*+f^,,(2f/,.+f/j (35) 

Notice that, although we have imposed planar symmetry, we have expressed the Ricci tensor in terms of 
derivatives of the conformal metric and of the conformal connection function just as is done for 
full 3-D numerical relativity 



C. von Neumann stability analysis about flat space 



As outlined for the scalar wave equation in Section II A, we discretize both the ADM equations (Eq 



and the CT equations (Eq. |3l|) using a 2-iteration iterative Crank-Nicholson scheme. The resulting finite 
difference equations, even though planar symmetry and a simplified form of the metric is assumed, are 
still too complicated to perform a von Neumann analysis by hand. The complication arises due to the 
recursive nature of the iterative Crank-Nicholson method; a 2-iteration procedure results in the source 
terms of Eqs. |2^ and ^ being computed 3 times in recursive succession. 

We have performed the von Neumann analysis of the ADM and CT equations in two independent ways. 
In the first way, we used the symbolical calculation computer program MATHEMATICA to explicitly 
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calculate the 2-iteration iterative Crank- Nicholson update, and then explicitly calculated the amplification 
matrix Eq. ^ in terms of the initial data variables. We then took these expressions, substituted the 
initial data about which we want to compute the von Neumann analysis, and calculated the eigenvalues 
to arbitrary precision inside MATHEMATICA. The second, independent way of performing the von 
Neumann analysis was to write an evolution code for the finite difference equations. We then input the 
initial data about which we want to compute the von Neumann analysis and computed the derivatives in 
the definition of the amplitude matrix Eq. |^ using finite differencing (this amounts to finite differencing 
the finite difference equations!). The package EISPACK was then used to compute the eigenvalues 
of the resulting amplification matrices. To obtain the highest accuracy possible, whenever calculating 
the amplification matrix Eq. ^ using finite differencing, we finite difference with multiple discretization 
parameters h, and use Richardson extrapolation to obtain values of the derivatives. 

Both methods were used, and found to produce identical results. All results reported in the paper were 
produced using both methods as described above. We emphasize the use of two independent methods 
not only to verify our results, but also for more practical reasons: one may eventually want to perform a 
von Neumann analysis for a full 3-D numerical relativity code. The analytic method using a symbolical 
manipulation package such as MATHEMATICA may not be feasible in the near future. It took well over 
100 hours on one node of an Origin 2000 running MATHEMATICA to perform the symbolical calculations 
needed for the von Neumann analysis of a plane symmetric code. It may be several orders of magnitude 
more expensive to analyze a full 3D evolution update. The finite difference method for computing the 
amplification matrix is much quicker, and it is reassuring that the finite difference method, used in 
conjunction with Richardson extrapolation, is accurate enough to reproduce the same detailed structures 
of the eigenvalues of the amplification matrices as doing the analytic calculation. 
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FIG. 2. The largest eigenvalue (ignoring eigenvalues that are identically 1) of the amplification matrix as a 
function of 6k = Az for the 2-iteration iterative Crank-Nicholson discretization of the ADM equations with 
flat initial data. Various values of the Courant factor a — At/ Ax are shown. 



In Figure |^ we plot the maximum of the modulus of the eigenvalues (neglecting eigenvalues that are 
exactly 1) of the amplification ma trix of the 2-iteration iterative Crank-Nicholson discretization scheme of 
the ADM equations from Section HI A using flat space as initial data, namely, gu = a = 1, and Ku = 0. 
We see that the spectral radius of the amplification matrix is less than or equal to 1 for a ~ At/ Ax < 1. 
Notice that for the Nyquist frequency mode {9k = kz Az = tt), when the Courant factor ct is 1, all 
eigenvalues of the amplification matrix have a modulus of exactly 1. For cr > 1, all Nyquist frequency 
modes are unstable (i.e. they all have amplification matrices with spectral radii > 1). 

In Figure |^ we plot the maximum of the modulus of the eigenvalues (neglecting eigenvalues that are 
exactly 1) of the amplification matrix of the 2-iteration iterative Crank-Nicholson discretization scheme of 
the CT equations from Section [II B, again using flat space as initial data. The resulting plot is identical 
to that of the ADM equations for ffat space. The stability properties of the CT equations about ffat 
space are exactly the same as the stability properties of the ADM equations about ffat space. 
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FIG. 3. The largest eigenvalue (ignoring eigenvalues that are identically 1) of the amplification matrix as a 
function of Ok = Az for the 2-iteration iterative Crank-Nicholson discretization of the CT equations with flat 
initial data. Various values of the Courant factor a — At/ Ax are shown. 



IV. VON NEUMANN ANALYSIS FOR NONLINEAR PLANE WAVES 

In this section, we study the stability properties of the ADM and CT equations about nonlinear plane 
waves. We require initial data that corresponds to nonlinear plane waves that satisfy the constraints and 
takes on the form of our simplified metric Eqs. ^ and We choose an exact plane wave solution first 
given by [p^ . The metric is assumed to take the form 

ds^ = -dt" + L^e^^dx" + e-2/5dy2) + dz'' , (36) 

where L = L{v), (3 = (3{v), and v — t — z. Given an arbitrary function /3(w), the Einstein equations 
reduce to the following ordinary differential equation for L(v): 

d'^L{v) fd/3{v)^^ 



dv^ n^j^^^^^' ^''^ 

In this paper, we take (3{v) to be given as 

f3{v) = ^e-^'^/i'^ (38) 

Eq. ^ is then solved with a 4th order Runge-Kutta solver. Initial data is obtained by setting t — (and 
thus, V ~ —2^), shown in Figure ^. We present results of a von Neumann analysis about the point z — 3.0 
with Az — 0.15. We have investigated different values of z and Az, and find the results presented to be 
generic for different values of z and Az, as long as we remain in the nonlinear regime, z < 8. 
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FIG. 4. Nonlinear plane wave initial data. Plotted are the metric functions Qxx and gyy as well as the extrinsic 
curvature functions K^x and Kyy. The metric function g^z as well as the lapse a are identically 1, and the 
extrinsic curvature component Kzz is identically 0. 

The main difference between the stability properties of discretizations of the CT and ADM systems 
is shown in Figures ^ and ^. These show, respectively, plots of the maximum modulus of the eigenvalues 
(ignoring eigenvalues that are exactly 1) of the amplification matrices for discretizations of the ADM and 
CT systems at 2 = 3.0. These plots show the range of the Courant factor 0.3 < a = Ai/Az < 0.9 and 
the range of modes O.Ttt < 9^ = Az < tt. Notice how, for these ranges of Courant factor a and mode 
wavenumber k^, all of the eigenvalues for the CT system are less than or equal to 1, where as for the 
ADM system, there is always at least one eigenvalue that has a modulus that is greater than 1. For these 
ranges of a and k^, we conclude that the CT system is stable, while the ADM system is unstable. 




k 

FIG. 5. The largest eigenvalue (ignoring eigenvalues that are exactly 1) of the amplification matrix as a 
function of 6k = kz Az for the Crank-Nicholson discretization of the ADM system about the nonlinear plane wave 
initial data at z = 0.3 with Az = 0.15. Different values of the Courant factor a = At/ Az are shown. 
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FIG. 6. The largest eigenvalue (ignoring eigenvalues that are exactly 1) of the amplification matrix as a 
function of Ok = Az for the Crank-Nicholson discretization of the CT system about the nonlinear plane wave 
initial data at z = 0.3 with Az = 0.15. Different values of the Courant factor a = At/ Az are shown. 



By looking at Courant factors of values a > 1, we see the typical Nyquist frequency instability, shown 
in Figures |^ and ^ for the ADM system and CT system, respectively. Notice in both Figures ^ and || 
that the largest modulus of the eigenvalues of the amplification matrices for long wavelength modes 
{dk = Az < O.Stt) are all greater than 1. In fact, this is the case for all values of Courant factor a. 
This is due to the fact that there is an exponentially growing gauge mode in the analytic solution to 
the analytic equations. Recall that we are using a different gauge choice than that given by the exact 
solution in Eqs. Using the gauge choice given in Eq. ^ and Eq. there exists an exponentially 

growing gauge mode in K^z- To take into account equations that admit exponentially growing solutions, 
the von Neumann condition, Eq. H, must be modified (see for details). In order that finite difference 
discretizations of equations that admit solutions that have exponentially growing modes remain stable, 
we must have 



\X^\ <l + 0{At). (39) 

To verify that the long wavelength {9k = kz Az < O.Stt) phenomena observed in Figures ^ and ||is simply 
due to the existence of (long wavelength) exponentially growing modes in the analytic solution to the 
analytic equations, we repeat the calculations leading to Figures ^ and ^, but decrease the discretization 
parameter Az by a factor of 2 (Az = 0.15/2 = 0.075). The resuhs for both the ADM and CT systems 
are similar, so we only show the results for the CT system in Figure ^. Notice that by decreasing Az 
by a factor of 2 and holding the Courant factor a constant, we also decrease At by a factor of 2. By 
comparing the long wavelength {6k = kz Az < O.Stt) sections of Figures ^ and we indeed see that the 
difference between the maximum modulus of the eigenvalues of the amplification matrices and 1 decreases 
by a factor of 2 when At is decreased by a factor of 2. Therefore, the discretizations of both the ADM 
and CT systems are stable for the long wavelength, exponentially growing gauge mode. Of course, the 
maximum eigenvalues for the high frequency sections of Figures ^ and ^ for ct > 1 do not approach 1 as 
0{At), which signifies a true von Neumann instability for cr > 1. 
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FIG. 7. The largest eigenvalue (ignoring eigenvalues that are exactly 1) of the amplification matrix as a 
function of = kz for the Crank-Nicholson discretization of the ADM system about the nonlinear plane wave 
initial data at « = 0.3 with = 0.15. Values for the Courant factor 0.9 < a < 1.05 are shown. 
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V. DISCUSSION AND CONCLUSIONS 



First, while the results presented in this paper are specific to the discretization method and initial 
data chosen here, one of the main points of this paper is to show that a von Neumann analysis can indeed 
be applied directly to discretizations of the Einstein equations. In the past, particular discretization 
methods were only analyzed with the von Neumann method through simple equations, such as the linear 
wave equation. Here we show that it is possible to carry out a von Neumann analysis on discretizations 
of equations as complicated as the Einstein equations. We would like to point out that a von Neumann 
analysis could also be used to test and/or formulate boundary conditions. The stability of outer boundary 
conditions, as well as the stability of inner boundary conditions for black hole evolutions, could be tested 
first through a von Neumann analysis instead of the traditional (and painful) method of coding up 
an implementation and looking for (and usually finding) numerical instabilities during the numerical 
evolution. 

In this paper, we have shown that the stability properties, as determined by a von Neumann stability 
analysis, of a common discretization (a 2-iteration iterative Crank-Nicholson scheme) of the ADM and 
CT systems about flat space are similar to the stability properties of the scalar wave equation. However, 
as we would like to emphasize again, the stability properties of a nonlinear finite difference update 
operator depend on the values of the discrete evolution variables. Therefore, it is not enough to study 
the stability of numerical relativity codes about flat space. In principle, one must verify the stability 
of a nonlinear finite difference update operator about every discrete state encountered during the entire 
discrete evolution process, as argued in Section ||. As a first step in this direction, we studied the 
stability properties of the ADM and CT systems about highly nonlinear plane waves by performing a 
von Neumann analysis for these scenarios. Several interesting features presented themselves. 

The main difference between the stability of the ADM and CT systems about the nonlinear wave 
solution is seen in Figures |^ and |6[ There, we see that, for a wide range of Courant factors a ~ At/ Az and 
wavenumbers (which includes the mode that is usually the most troublesome in numerical relativity, 
the Nyquist frequency mode whose wavelength is 2 Az), the CT system is stable (i.e., the amplification 
matrix has a spectral radius of 1 or less) whereas the ADM system is unstable (i.e., the amplification 
matrix has a spectral radius greater than 1). Note that the ADM system, in this region of Courant factor 
a and wavenumber k^, has an amplification matrix whose spectral radius is approximately 1 + 10~^. This 
is a very small departure from unity. Thus, instabilities arising from these modes could take a long time 
to develop during numerical evolutions. For exam ple, t he 2-iteration iterative Crank-Nicholson update 
operator for the scalar wave equation from Section II A has a spectral radius of 1 -I- 10^^ for a Courant 
factor a ~ At/Az = 2.000004, and this update operator can be used to evolve initial data sets many 
wavelengths before the Nyquist frequency instability sets in. 

Again, it must be pointed out that these results are specific to plane wave spacetimes, simplified 
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(diagonal) forms of the metric, choice of gauge, choice of initial data used here, and the 2-iteration 
iterative Crank-Nicholson scheme. It remains to be seen whether or not these results are generic to other, 
more general discretizations of the Einstein equations and for more general data, such as black hole and 
neutron star discrete evolutions. The one thing that is certain is that the von Neumann analysis can be 
used as a diagnostic tool for determining the stability of discretizations of nonlinear differential equations 
as complicated as the Einstein equations. 
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